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We develop a finite temperature Hartree theory for the trapped dipolar Bose gas. We use this theory to study 
thermal effects on the mechanical stability of the system and density oscillating condensate states. We present 
results for the stability phase diagram as a function of temperature and aspect ratio. In oblate traps above the 
critical temperature for condensation we find that the Hartree theory predicts significant stability enhancement 
over the semiclassical result. Below the critical temperature we find that thermal effects are well described by 
accounting for the thermal depletion of the condensate. Our results also show that density oscillating condensate 
states occur over a range of interaction strengths that broadens with increasing temperature. 
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I. INTRODUCTION 

A significant new area of interest in ultra-cold atomic gases 
is the study of systems in which the particles interact via 
a dipole-dipole interaction (DDI) [1]. This interest is be- 
ing driven by a broad range of proposed applications from 
condensed matter physics to quantum information, e.g. see 
[2]. Experimental progress in the quantum degenerate regime 
has been driven by seminal work with 52 Cr [3], which was 
Bose condensed in 2005, and more recently the realization of 
Bose-Einstein condensates of 164 Dy [4] and 168 Er [5]. Po- 
lar molecules, which have DDIs several orders of magnitude 
larger than those of the atomic gases, have already been pro- 
duced in their ground rovibrational state [6, 7], and steady 
progress is being made towards cooling these into the degen- 
erate regime. We also note the recent achievement of a degen- 
erate Fermi gas of 161 Dy [8]. 

The DDI is long-ranged and anisotropic with both attractive 
and repulsive components. Therefore, an important consider- 
ation is under what conditions the system is mechanically sta- 
ble from collapse to a high density state. Theoretical studies 
on zero temperature dipolar condensates reveal a rich stabil- 
ity diagram where, due to the DDI anisotropy, the stability is 
strongly dependent on the geometry of the trapping potential 
and the properties of the short ranged (contact) interactions 
[9-12]. Another interesting theoretical observation is that for 
appropriate parameters (near instability) the condensate mode 
exhibits spatial oscillations and has a density maximum away 
from the minimum of the trapping potential [11-14]. How- 
ever, evidence for this density oscillating state has yet to be 
observed in experiment. 

In this work we study the properties of a trapped dipolar 
Bose gas at finite temperature - a regime largely unexplored 
in theory and experiments. In previous work [15] we stud- 
ied the stability of a normal Bose gas (i.e. above T c ) using a 
self-consistent semiclassical approximation. In this work we 
extend this study to below T c and to include quantum pressure 
(i.e. beyond-semiclassical effects) by numerically solving for 
the condensate and its excitations. Using this theory we study 
the crossover from the high temperature (above T c ) to zero 
temperature (pure condensate) stability. Our results reveal that 
beyond semiclassical effects play a significant role above T c 
in oblate geometry traps and enhance the stability region, and 



that the double instability phase diagram in this trap geome- 
try (predicted in [15]) remains prominent. We also study the 
behavior of the emergent biconcave condensate (density oscil- 
lating ground state) in the finite temperature regime, and find 
that thermal effects enhance the density oscillation and en- 
large the parameter regime over which this type of state exists. 
We demonstrate that the below T c temperature dependence of 
the stability boundary is well-characterized by a simple model 
that accounts for the thermal depletion of the condensate. 



n. FORMALISM AND NUMERICAL IMPLEMENTATION 

A. Formalism 

Here we consider a set of particles of mass M confined in 
a cylindrically symmetric harmonic potential 

C/ tr (x) - \M[u%x 2 + y 2 ) + lo 2 z z\ (1) 

of aspect ratio A = oj z /ui p , with z and p representing the axial 
and radial directions, respectively. We take the particles to 
have dipole moments polarized along the z axis by an external 
field, such that the DDI potential between particles is 
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where Cdd = Mo A* 2 (=<^ 2 / e o) is the magnetic (electric) dipole- 
dipole interaction strength and 8 is the angle between the z di- 
rection and the relative separation of the dipoles (r = x' — x). 
It is easy to extend our calculations to include local (con- 
tact) interactions, however here we focus on the case of pure 
dipole-dipole interactions, as has been realized in experiments 
by use of a Feshbach resonance (e.g. see [16]). 

The Hartree formalism we employ (see Appendix A for a 
discussion of the relation to Hartree-Fock theory and relevant 
terms neglected) involves solving for the system modes using 
the non-local equation 
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where 



y cff (x) = C/ tr (x) + Jdx' f/ dd (x - x')n(x'), (4) 
n(x) = ^AT,Mx)| 2 , 



(5) 



are the effective potential and total density, respectively, with 
Nj = [e^' e '~w — 1] 1 the equilibrium (Bose-Einstein) oc- 
cupation of the mode, /? = the inverse temperature, 
and /i the chemical potential. Equations (3)-(5) are solved 
self-consistently while the chemical potential is adjusted to 
ensure that the desired total number N = j d 3 xn(x) is ob- 
tained. Below the critical temperature T c a condensate forms 
in the lowest mode un(x) with No ~ N, but the theory, as 
written in Eqs. (3)-(5), requires no additional adjustment to 
account for the condensate (due to our neglect of exchange) 
and smoothly transitions across T c . We emphasize that our 
motivation for using this theory is that it includes the domi- 
nant direct interactions and the full discrete character of the 
low energy modes, yet is more computationally efficient than 
Bogoliubov-based approaches. This enables us to study chal- 
lenging problems that have not been explored, in particular fi- 
nite temperature mechanical stability, in which obtaining con- 
vergent self-consistent solutions is demanding and time con- 
suming. Our numerical approach builds on various develop- 
ments (particularly those described in [17]) and includes a 
number of features to aid calculations in the finite tempera- 
ture regime where interaction effects dominate (see Appendix 
B for details). 

The neglect of dipole exchange is consistent with other 
work on finite temperature bosons [18] and zero temperature 
studies of fermion stability [19]. We would like to note that 
there is some justification for this approximation. Studies on a 
normal trapped dipolar Fermi gas suggest that exchange inter- 
actions will quantitatively, but not qualitatively, affect stabil- 
ity [19, 20]. Indeed, the thermodynamic study of that system 
presented in [21] found that exchange interactions are typi- 
cally less important than direct interactions except for traps 
that are close to being isotropic. Similarly, Ticknor studied 
the quasi-two-dimensional Bose gas using the Hartree-Fock- 
Bogoliubov-Popov (HFBP) meanfield theory [22] and found 
that exchange terms were generally less important than direct 
terms. 



III. RESULTS 
A. Comparison to previous calculations 

To benchmark our Hartree calculations we perform a quan- 
titative comparison to the HFBP calculations that Ronen et 
al. [18] performed for the three-dimensional trapped Bose gas 
at finite temperature. In Sees. Ill A 1 and III A 2 we make this 
comparison for two different sets of results from [18]. 

We note that those HFBP calculations excluded thermal 
exchange interactions, although they did include condensate 
exchange interactions (exchange interaction of condensate 



atoms on the thermal excitations) [23]. We extended our 
Hartree algorithm to include condensate exchange but found 
it made negligible difference to the predictions and do not in- 
clude results with this term here. 



1. Condensate Fraction 

The results of the first comparison we perform are presented 
in Fig. 1(a). There we compare the condensate fraction, as a 
function of temperature, for a system with A = 7. We ob- 
serve that the Hartree and HFBP theories predict an appre- 
ciably lower condensate fraction than the ideal case, and are 
in very good agreement with each other over the full tem- 
perature range considered. The low energy excitations of a 
Bose-Einstein condensate are quasi-particles, which are ac- 
curately described by Bogoliubov theory (such as the HFBP 
theory), however the thermodynamic properties of the system 
are dominated by the single particle modes (e.g. see [24]). A 
comparison of the Bogoliubov and Hartree-Fock spectra of a 
T = dipolar Bose-Einstein condensate (BEC) was made 
in [17]. That comparison revealed that the spectra were al- 
most identical, except for low energy modes with low values 
of angular momentum, where small differences in the mode 
frequencies were observed. 




FIG. 1. (a) Condensate fraction and (b) density oscillation contrast 
(see text) for a dipolar BEC in a A = 7 pancake trap. Hartree re- 
sults (pluses), HFBP results (solid lines), ideal gas result (dashed 
line). HFBP data corresponds to results shown in Figs. 5 and 6 of 
Ref. [18], Other parameters: {uj p ,uj z } = 2tv x {100, 700} s"\ 
N — 16.3 x 10 3 52 Cr atoms with contact interactions tuned to 
zero. T° = y N/((3)huj/kB is the ideal condensation tempera- 
ture, where lj — ^/ujpUj z and C( Q ) is the Riemann zeta function 
withC(3) » 1.202. 



2. Density Oscillating Ground States 

An interesting feature of dipolar condensates is the occur- 
rence of ground states with density oscillation features, where 
the condensate density has a local minimum at trap center. For 
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a cylindrically symmetric trap these states are biconcave (red 
blood cell shaped - surfaces of constant density are shown in 
Fig. 6) first predicted for T = condensates in Ref. [11]. 
In the purely dipolar case such biconcave states occur under 
certain conditions of trap and dipole parameters, but notably 
only for A > 6 and for dipole strengths close to instability. 
In [18] the HFBP technique was used to assess the effect of 
temperature on the density oscillating states. This was char- 
acterized by the contrast, a measure of the magnitude of the 
density oscillation, defined as 
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where n(0) is the density at trap center and n max is the maxi- 
mum density of the system. 

In Fig. 1(b) we compare our Hartree and HFBP theories for 
the contrast. This comparison reveals some small residual dif- 
ferences between the theories, however the results are in rea- 
sonable agreement and both predict that the contrast goes to 
zero (i.e. the condensate returns to having maximum density 
at trap center) at T m 0.65T C ° . 



B. Mechanical stability 

Our first application of the Hartree theory is to study the 
finite temperature mechanical stability of a trapped dipolar 
Bose gas. To do this we construct a phase diagram for the 
range of dipole strengths for which the gas is stable for a num- 
ber of different trap geometries. Such stability properties, and 
the dependence on interactions and trap geometry, have been 
measured accurately in the dipolar system in the zero temper- 
ature limit (e.g. see [16]). We note theoretical studies [25-27] 
showing the important role of temperature on the observed 
stability of 7 Li condensates [28], which have an attractive con- 
tact interaction. 



1. Locating the stability boundary 
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FIG. 2. Locating instability (upper subplots): (a) The total num- 
ber of atoms of the self-consistent Hartree solution versus chem- 
ical potential for A = 1/8, fcsT = AOhui and Cdd = 7 x 
10 -4 fojah o (dashed, case A), 1.5 x 10" 4 fiwaL (solid, case B), 
with aho = y/K/Mu). Each line terminates at the point of in- 
stability and occurs at the respective critical number iVcrit. (b) 
Same results as in (a) but plotted against eq — /Li. The dotted line 
represents the target number, in this case N = 2 x 10 . Den- 
sity Profiles (lower subplots): Solid (dashed) line represents the ra- 
dial (axial) density n, higher curves are near the stability bound- 
ary. A = 1, N = 2 x 10 5 . (c) T/T c ° = 0.82 (N /N « 0.43) 
and Cdd/47rC = {3.65 x 10" 4 (gray), 1.83 x 10~ 4 (black)}, 
(d) T/T c ° = 1.27 and C di /47vC = {2.91 (gray), 1.22 (black)}. 
We have introduced the interaction strength unit Co = Hua^ / yN 
which is convenient for cases where iV is fixed, and allows our sub- 
sequent results to be directly compared to those in [15]. 



We consider a trapped sample of fixed mean number N and 
wish to determine the values of the dipole interaction param- 
eter for which the system is mechanically stable as a function 
of temperature. In doing so we construct a phase diagram 
in {Cdd, T}-space that indicates the stable region. In prac- 
tice we locate the stability boundary (i.e. a curve) that sep- 
arates the stable and unstable regions. Our procedure to ob- 
tain this boundary involves a computationally intensive search 
through parameter space to find the self-consistent solutions 
on the verge of instability. Determining the stability boundary 
for fixed mean number N complicates this process: since we 
work in the grand-canonical ensemble where the proper vari- 
ables are {/i, Cdd,? 1 }, an additional iterative search over the 
parameter /i is required to fix N to the desired target number. 

In Fig. 2 we provide some examples to illustrate how we 
identify the value of the DDI at the stability boundary for a 
gas with (target number) N — 2 x 10 5 atoms at a particular 



temperature. To do this we show the dependence of total atom 
number on /i for two different values of the DDI [Fig. 2(a)]. 
For both curves the total number increases as we move along 
these curves until some maximum value -/V cr i t is reached at 
which the system becomes unstable. The non-monotonic be- 
havior of these curves arises because the ground state energy 
e changes as the number of atoms increases, and hence the 
role of DDIs increases. For this reason we also show the same 
two cases, but as a function of eo — M> m Fig- 2(b). 

The sharp cusps in Figs. 2(a) and (b) correspond to the point 
where the system condenses [i.e. where eo — fi ~ 0]. The 
dependence of eo on 7V n is strongly dependent on the trap ge- 
ometry, and for the cases we consider here with A = 1/8, eo 
decreases with increasing Nq. This is because the head-to- 
tail character, in the cigar geometry, emphasizes the attractive 
part of the DDI so that as the condensate number increases, eo 
(« n) decreases. 
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For case A in Fig. 2(a) the number at which collapse occurs 
is less than the target number, thus we conclude that the DDI 
used in this calculation lies within the unstable region for the 
system (i.e. no stable solution can be found for N atoms with 
this value of DDI). In contrast, for case B in Fig. 2(a) N cr a > 
TV and thus the value of DDI is in the stable region. To locate 
the stability boundary points we need to trace out these curves 
for various Cdd values until we find N a it = N to within our 
numerical tolerance (this has to be done for each value of T). 
This process is painstaking and can take several days to find a 
single point on the stability boundary. 

We identify the self-consistent Hartree solutions as being 
unstable when they become grid-dependent. This means that 
as the distance between grid points tends to zero, the radial 
width of the cloud contracts and the chemical potential tends 
to negative infinity. Precisely locating the instability point is 
a stringent numerical task and requires careful convergence 
tests. For condensates with contact interactions this type of 
numerical instability analysis was applied in Refs. [25-27] 
(also see Ref. [15]). In Fig. 2(a) the instability point occurs 
at the end of the upper horizontal plateau in the N versus fi 
curves (compare to Fig. 1 of [25]). We show examples of the 
spatial density profiles for a spherical trap in Figs. 2(c) and 
(d). The system considered in Fig. 2(c) is condensed, while 
that considered in Fig. 2(d) is above the critical temperature. 
For both cases a result is shown that is well inside the stable 
region (black curves) and near the stability boundary (gray 
curves). Despite a large difference in the density scales of 
the two regimes they both exhibit a similar sharpening of the 
density profile near instability. 

An additional consideration emerges for stability calcula- 
tions below T c in regimes where the condensate is in a density 
oscillating state. Here the first mode to go soft (and then de- 
velop imaginary parts) as the stability boundary is reached is a 
m 7^ quasi-particle mode [29], where m is the angular mo- 
mentum projection quantum number (so called angular roton 
mode [11]). This instability is not revealed in the Hartree exci- 
tations, and as we solve for the condensate in the m = space 
(see Appendix B), the condensate does not exhibit numerical 
instability. Thus in cases where the condensate exhibits a den- 
sity oscillating state we perform a Bogoliubov analysis of the 
condensate mode (within the effective potential of the self- 
consistent Hartree solution) to determine if any ra/0 modes 
have become unstable [30]. 

2. Stability above T c 

In Fig. 3 we show our results for the stability of the normal 
phase. In previous work we examined this regime using a 
semiclassical Hartree approach in which the density is 

n(x) = A^ 3 C3/ 2 (e^- v -Ml) , ( 7 ) 

where V c g (x) is the effective potential calculated using n(x) 
[see Eq. (4)], Ca( z ) = YjjLi ^ /j a i s me B° se function, and 
AdB = h/y / 2irAlkBT. The semiclassical results are shown 
as solid lines in Fig. 3. 




FIG. 3. (Color online) Stability regions in DDI-temperature space. 
Shaded regions indicate stability for each geometry, from top to bot- 
tom A = {8, 4, 2, 1, 1/2, 1/8), the geometric mean trap frequency is 
fixed and N = 2 x 10 s . Actual data points represented by symbols 
while the shading of the stable regions interpolates to guide the eye, 
the semiclassical model is given by the solid curves. Error bars rep- 
resent the 1 a spread in the convergence test (see Appendix B 3 for 
more details). 

We observe that as a general trend the stability region grows 
with increasing A. The strong geometry dependence of these 
results arises from the anisotropy of the dipole interaction: In 
oblate geometries (A > 1) the dipoles are predominantly side- 
by-side and interact repulsively (stabilizing), whereas in pro- 
late geometries (A < 1) the attractive (destabilizing) head-to- 
tail interaction of the dipoles dominates (a similar geometry 
dependence is observed for the stability of T = dipolar con- 
densates [11, 12]). 

A primary concern is the nature of beyond semiclassi- 
cal effects, i.e. what differences emerge from our diagonal- 
ized Hartree theory over the semiclassical formulation. Most 
prominently in the results of Fig. 3 we observe that while 
the Hartree and semiclassical stability boundaries are in good 
agreement for prolate geometries, in oblate traps the Hartree 
results are significantly more stable. This difference between 
the boundaries predicted by the two theories increases with 
increasing A. This observation is surprising because our cal- 
culation is for a rather large number of atoms (iV = 2 x 10 5 ), 
where the semiclassical approximation would normally be ex- 
pected to furnish an accurate description of the above T c be- 
havior. 

We attribute this failure of the semiclassical theory to its 
inappropriate treatment of the interactions between the low 
energy modes [31]. The nature of the DDI, when tightly con- 
fined along the polarization direction, has been extensively 
studied in application to pure BECs [14, 32], where it has 
been shown that it confers additional stability on the system, 
as verified in recent experiments [33]. This arises from a con- 
finement induced momentum dependence of the interaction: 
the interaction is repulsive (stabilizing) for low momentum 
interactions, but decays to being attractive with a character- 
istic wavevector k ~ 1 ja z set by the z confinement length 
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a z = y/h/Muj z . Notably these features of the confined inter- 
action mediate BEC instability through the softening of radi- 
ally excited modes with a wavelength ~ a z [32, 34—37]. 

It is not clear that these confinement effects will be applica- 
ble at a modestly oblate trap with A = 8, however numerical 
studies have revealed that quasi-particle modes with a wave- 
length ~ a z soften in a BEC with A = 7 [11]. Within the 
limited range of results we have for A > 1 we see evidence 
consistent with confinement induced effects playing an impor- 
tant role in the above T c Hartree calculations. Notably, that 
the relative difference between the stability boundaries of the 
Hartree and semiclassical calculations scale with l/a 2 z . Also, 
when the system is unstable, during the self-consistency it- 
erations (prior to collapse) strong radial density fluctuations 
develop in the system 

A key prediction from our semiclassical study [15] is a dou- 
ble instability feature in oblate trapping geometries arising 
from the interplay of thermal gas saturation and the anisotropy 
of the DDI. Our Hartree calculations in this oblate regime, 
despite shifting the stability boundary from the semiclassical 
prediction by a considerable amount, reveal that the double 
instability feature is robust to beyond-semiclassical effects. 

A prominent feature of the semiclassical calculation is that 
the stability curves for the purely dipolar gas terminate at the 
critical point with Cdd = (i.e. predicting that without con- 
tact interactions only an ideal gas is stable below T c ). This oc- 
curs because the local compressibility at trap center diverges 
at the critical point and the gas is unstable to any attractive in- 
teraction (see [15]). In the beyond-semiclassical calculations 
the trap provides a finite momentum cutoff that prevents the 
divergence of compressibility, and thus the system has a fi- 
nite residual stability at and below T c (which we consider in 
Sec. Ill B 3). 



3. Stability below T c 

In Fig. 4 we consider the stability below T c where the semi- 
classical model does not apply. These results are identical to 
those shown in Fig. 3, but the below T c details are revealed us- 
ing a logarithmic vertical axis. Compared to the above T c gas 
the condensate is rather fragile, with the critical DDI strength 
defining the stability boundary decreasing by ^3 to 4 orders 
of magnitude. 

In the zero temperature limit our results agree with previ- 
ous calculations based on solving the Gross-Pitaevskii equa- 
tion [11]. This agreement is expected as the two theories are 
identical when the excited modes have vanishing population. 
For a pure condensate, the critical DDI strength depends on 
the condensate number and trap geometry according to [11] 



C d * d = ^P, (T = 0) 



(8) 



with F( A) a rather interesting function of trap geometry alone, 
as characterized in Fig. 1 of [1 1] [41]. More generally, beyond 
the case of pure DDIs, F also depends on the contact interac- 
tion strength, e.g. see [12, 37]. 




FIG. 4. (Color online) Stability boundary focusing on the below T c 
behavior (line styles as in Fig. 3). For reference the ideal finite size 
adjusted critical temperature T Ci fs [38] for two geometries (T^Js 
and T^p%) are indicated by short vertical lines. The effect of DDIs 
on T c was calculated perturbatively in [39, 40], however our results 
are far outside the perturbative regime. 
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FIG. 5. (Color online) Stability boundary scaling. The stability 
boundary results (symbols) have been taken from Fig. 4 for A = 
{8,1,1/8} (top to bottom). Dashed line prediction is based on a non- 
interacting iVo scaling (see text) and the solid line uses the iVo calcu- 
lated from the Hartree solutions. 



As temperature increases, but focusing on T < T®,we ob- 
serve in Fig. 4 that the stability boundary increases signifi- 
cantly. This occurs because as the temperature increases the 
condensate is thermally depleted. Indeed, by simply account- 
ing for the thermal depletion we can immediately extend result 
(8) to predict the critical value of the DDI at finite temperature 



CUT) 



F(X) 
N (T) 



cS,(o) 



N 
JVo(T)' 



(9) 



where the last expression is obtained using N (T = 0) = N. 
Equation (9) predicts that the stability at finite temperature 
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increases inversely proportional to the condensate occupation 
and, as shown in Fig. 5, provides a good description of the sta- 
bility predictions from the full Hartree calculations. In these 
comparisons we have used two models for the condensate oc- 
cupation: (i) the non-interacting prediction 

iV NI (T) = N[l - (T/T°) 3 }, (10) 

and (ii) the value of No (T) obtained from the Hartree calcu- 
lations. Equation (9) using A/g rl (T) provides a good predic- 
tion for Hartree stability curves with A = 1/8 and 1. For the 
oblate system (A = 8) agreement is not as good as is apparent 
in Fig. 4 for T > 0.5T°. In this case the values of Cdd at the 
stability boundary are much higher than for the other values 
of A, and thus interaction effects more significantly affect the 
condensate. However, much better agreement is obtained if 
we take No (T) from the Hartree solution. 

We note that the simple model (9) does not account for any 
other effects of the thermal cloud on the condensate [e.g. ther- 
mal back action through modifications of V^g (x)]. Thus, the 
level agreement of this simple model with the full Hartree re- 
sults suggest that these additional effects are not significant in 
the regimes we have studied. 



C. Thermal effects on biconcavity 

As our final application we consider thermal effects on the 
biconcave density profiles which occur in oblate geometries. 
To date, the only study of temperature effects of these states 
was inRef. [18] [which we reproduce in Fig. 1(b)]. That study 
considered a single line (at fixed Cdd and N and varying T) 
through the phase diagram, and showed that biconcavity per- 
sisted at small finite temperatures (T < 0.25X)?), but then was 
rapidly washed out as temperature increased further. 

Using our Hartree theory we provide a broad characteriza- 
tion of the thermal effects on biconcavity. We focus on the 
case A = 8, which supports a biconcave condensate at T = 0. 
In Fig. 6(a) we present contours of biconcave contrast [as de- 
fined in Eq. (6)] over the entire range of parameters where 
this state is stable. These results show that biconcavity is not 
destroyed as temperature increases. Instead the parameter re- 
gion over which biconcavity occurs grows, with large bicon- 
cave contrasts emerging at higher temperature. The general 
trends seen can be understood by considering the thermal de- 
pletion of the condensate, using similar arguments to those 
made to obtain Eq. (9): as the temperature increases the value 
of Cdd required for the condensate to exhibit a biconcave den- 
sity profile should increase in a manner that is approximately 
inversely proportional to the condensate occupation. Thus, the 
washing out observed in [18] [our Fig. 1(b)] arises because 
they considered Cdd fixed. Thermal depletion of the conden- 
sate is not sufficient to explain all aspects observed in our re- 
sults, e.g. the deepening of the biconcave contrast that devel- 
ops at higher temperatures in Fig. 6(a). This arises from ad- 
ditional effects of the thermal interaction with the condensate, 
e.g. small changes in the aspect ratio of the effective poten- 
tial that the condensate experiences can significantly change 
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FIG. 6. (Color online) Biconcave characteristics for A = 8 and 
N = 2 x 10° at finite temperature, (a) Stability diagram with bi- 
concave contrast contours (0, 0.05, 0.1, 0.15, 0.2, 0.25} (bottom 
to top) added. The solid curves are interpolations between the cal- 
culated contour points. The white dotted line marks where we ter- 
minate the contours due to the condensate fraction becoming negli- 
gibly small. Triangles indicate the stability boundary from Fig. 4. 
Inset: Magnification of the high temperature region, (b) Radial 
densities for phase space points marked by A and B in (a). A: 
T/T° = 0.0910, Cdd/47rCo = 0.00268 and N /N = 1.00 (ther- 
mal depletion < 1%). B: T/T c ° = 0.910, C dd /4.7rC = 0.0291 and 
No/N = 0.0716. Solid (dashed) lines represent the total (conden- 
sate) density. Insets: corresponding surface contours at 67% of the 
peak density. 



the contrast (c.f. the strong dependence of biconcavity on trap 
aspect ratio near A = 8 in Fig. 1 of Ref. [11]). 

In Fig. 6(b) we show two examples of the biconcave density 
profiles at different temperatures. Case B displays the very 
pronounced biconcavity for a system at T m 0.9T®, where 
the condensate fraction is N /N s» 0.07. 



IV. CONCLUSIONS 

In this paper we have developed a Hartree theory for a 
trapped dipolar Bose gas that can be applied to make predic- 
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tions above and below the condensation temperature. We have 
used this theory to quantify the role of thermal fluctuations on 
the mechanical stability of the cloud, and present results for 
the stability phase diagram as a function of temperature and 
aspect ratio. Above T c our theory predicts significant correc- 
tions to the stability boundary over the equivalent semiclas- 
sical theory. Most notably, the semiclassical theory underes- 
timates the size of the stability region for oblate geometries. 
Below T c (but at finite temperature) we find that the stability 
boundary is well described by the zero temperature result after 
scaling for the thermally depleted condensate. 

We have also studied the role of thermal fluctuations on 
biconcave condensate states. Our results show that as tem- 
perature increases, and the condensate thermally depletes, the 
range of interaction parameters in which these kinds of states 
can be found increases. Furthermore, a large thermal cloud 
may actually enhance the biconcave contrast making direct 
imaging of an in situ density oscillating state more feasible, 
see Fig. 6(b). 

An important step for future theoretical studies in the fi- 
nite temperature regime is to fully include thermal exchange 
effects. Because a large number of modes are important 
for temperatures of the order of T c , full Hartree-Fock cal- 
culations will probably not be feasible. It is possible to in- 
clude exchange interactions semiclassically (c.f. Fermi stud- 
ies [20, 21]), although our work here has shown beyond semi- 
classical effects are important even above T c . An interesting 
possibility is the extension of classical field methods to ther- 
mal dipolar gases (e.g. [42—45]) which may provide a compre- 
hensive description for temperatures around the condensation 
transition. 
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Appendix A: Hartree and Hartree-Fock theory for dipolar Bose 

gases 

In this appendix we describe the full Hartree-Fock the- 
ory for the dipolar Bose gas and discuss the reduction to the 
Hartree form we employ here. We then introduce the semi- 
classical Hartree approach we use to calculate high energy 
modes, which are insensitive to beyond-semiclassical effects. 



a. Above T c 

The Hartree-Fock equation for the modes of an uncon- 
densed dipolar Bose gas is 

e,«,(x) = H 0Uj ( X ) + / rf V t / dd (x'-x W x> 3 (x) 

V ' 

Hartree/Direct interaction term 

+ |d 3 x'[/ dd (x'-x)G(x,x> J (x'), (Al) 

S v ' 

Fock/Exchange interaction term 

where H is the single particle Hamiltonian and 

G(x,x') = ^ A/^*(x>,(x), (A2) 



n(x) = G(x,x), 



(A3) 



are the first order coherence function and total density, respec- 
tively. 



b. Below T c 

Below T c an appreciable number of atoms condense into 
the lowest energy single particle mode, which we denote as 
the condensate mode u (x) with respective energy e an d oc- 
cupation A ~ O(N). In this regime the Hartree-Fock equa- 
tions take the form 

e j u j (x) = H u j (x) + J d 3 x'C/ dd (x' - x)n(x>j(x) (A4) 
s ' 

condensate + thermal direct 

+ J rf 3 x' ( 7 dd (x'-x)G(x,x')^(x') 

s v y 

thermal exchange 

+ Q y dV [7 dd (x' - x)G (x',x)QK (x')}}, 



condensate exchange 



where 



G (x',x) = JV uS(x / )u (x) ) (A5) 
G(x',x) = AH(x'K(x), (A6) 



n(x) = G (x, x) + G(x, x), 



(A7) 



1. Hartree-Fock equations 

The Hartree-Fock theory for a Bose gas is well-established 
[46], particularly for the case of contact interactions (e.g. see 
Refs. [25, 26]). Here we present this theory for a system inter- 
acting with a DDI and consider the cases of above and below 
T c separately. 



are the condensate and thermal first order coherence functions, 
and the total density, respectively. We have also introduced the 
projector 

Q{/(x)} = J d 3 y [<5(x - y) - uo(x)«o(y)] /(y), (A8) 

to remove components of /(x) parallel to the condensate 
mode. The projection operator in Eq. (A4) acts to ensure that 



atoms within the condensate do not undergo an exchange in- 
teraction with themselves. In particular, when acting on the 
condensate mode Eq. (A4) reduces to the expected general- 
ized Gross-Pitaevskii equation 



(a) 



(b) 



e w (x) 



H + J d 3 x'£/ dd (x'-x)n(x') 



uo(x) (A9) 



tfVt/ dd (x'-x)G(x',x)u (x'), 



which has direct interactions with condensate and thermal 
atoms, but only thermal exchange. The projector also ensures 
that the modes {uj(x)} form an orthonormal set (e.g. see 
[47,48]) [49]. 



2. Reduction to Hartree theory 

The full numerical solution of Eq. (A4) [or even Eq. (Al)] 
is extremely challenging, because the thermal exchange re- 
quires G(x', x) to be calculated (or at least applied) for each 
self-consistent iteration. This limits the theory to applications 
involving a small number of modes and away from regimes 
where mechanical stability can be studied. 

The Hartree theory we use is obtained by neglecting 
both condensate and thermal exchange terms [as labeled in 
Eqs. (A4)], yielding 



w( x ) 



Hq + J d 3 x' £7 dd (x' - x)n(x') 



Uj(x), 



(A10) 



[i.e. Eq. 3]. The absence of exchange terms means that a pro- 
jector is no longer needed. 

The properties of the Hartree, Hartree-Fock and other theo- 
ries for the homogeneous Bose condensed gas, including con- 
servation laws, are extensively discussed in Sec. VI Ref. [50] 
(also see [51] for a discussion of the HFB and HFBP theories 
of the inhomogeneous system). We note that in the uniform 
purely dipolar gas, the Hartree term is zero and the DDIs af- 
fect the system only though the Fock term. However, in the 
trapped system the Hartree term is often dominant, particu- 
larly when the harmonic trap is appreciably anisotropic (i.e. 
A > 1 or A < 1). 



Appendix B: Description of Hartree algorithm 

In this section we discuss our implementation of the Hartree 
model as a numerical algorithm. 





FIG. 7. (Color online) Schematic of our approach to solving dipo- 
lar Hartree equations, (a) The low energy (discrete) modes are ex- 
plicitly solved for by diagonalizing the Hartree equations, whereas 
higher energy modes are treated using a semiclassical approxima- 
tion, (b) Schematic of the range of the spatial grids used to solve 
for the discrete modes and the large grid used to solve for the high 
energy semiclassical theory. 



modes. Instead we employ a hybrid method and diagonalize 
for the lowest energy modes, up to some energy e cut , and then 
calculate the remainder within the semiclassical approach [see 
Fig. 7(a)]. 

The semiclassical approach can be obtained by making the 
replacement V —> ik in Eq. (3), where k is a wavevector. 
This transforms the Hartree equation to an algebraic equation 
in which the energy is given in (x, k) -phase-space as 

e(x,k)= ^ + T/ eff(x ), (Bl) 

where the effective potential was defined in Eq. (4). The semi- 
classical portion of the system is described by a (Wigner) dis- 
tribution function W(x,k) = { e /3 [ £(x ' k)_ ' 1 ] - l} -1 . For con- 
sistency the semiclassical description can only be applied to 
regions of phase-space where e(x, k) > e cut to avoid dou- 
ble counting of modes. From this we obtain the semiclassical 
region density 

r (731, 

n sc (x)= / M/(x,k), (B2) 



1 



where ( 3/2 (z,y) = (2/Vtt) f y 
complete Bose function and 



C3/2 (e^- v ^\f3K min (x) 



(B3) 



\e l /z - ly-y/idt is the in- 



K min (x) = max{e cut - V^x), 0}. 



(B4) 



1. Semi-classical treatment of high energy modes 

The Hartree equation (3) is cylindrically symmetric and can 
be solved using a set of two-dimensional grids. However, at 
finite temperature typically > 10 5 modes are thermally acces- 
sible in the regimes we study, and a full self-consistent calcu- 
lation is not feasible in terms of the discrete (i.e. diagonalized) 



2. Summary of algorithm and numerical considerations 

All the excitations up to a given energy e cut are solved for 
using the Arnoldi algorithm. The associated discrete mode 
density is constructed 

n d (x)= ]T iV>i(x)| 2 , (B5) 

«j<«cut 
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The semiclassical and the total densities are then evaluated 

n(x) = n d (x) + n sc (x). (B6) 

We use fixed point iteration of these steps to ensure self con- 
sistency. To avoid oscillations only a small amount of the new 
prediction for the total density (n now ) is mixed in with the 
existing prediction (n id), i- e - 

n(x) -» A m i x n ncw (x) + ( )n id(x), (B7) 

where A m ; x is the mixing parameter. Upon obtaining a self- 
consistent solution external parameters are adjusted to tune 
the solutions to a desired macrostate (e.g. an outer loop of p, 
being adjusted to obtain the correct total number TV). 

We briefly mention a number of aspects of our algorithm: 

1. We make use of Fourier-Hankel techniques [17] to uti- 
lize the cylindrical symmetry and reduce the eigenvalue 
problem to being two-dimensional. The Fourier-Hankel 
approach is useful because it allows accurate Fourier 
transforms to simplify the evaluation of the convolu- 
tion required to construct the direct dipolar interaction. 
However, the radial Hankel transform requires a dif- 
ferent radial grid for each angular momentum projec- 
tion quantum number m, thus the problem requires a 
set of two-dimensional grids (we typically diagonalize 
modes with m up to 10, i.e. requiring 11 grids - gener- 
ally even more in oblate geometries). This requires ex- 
tensive transformation of quantities [e.g. n(x)] between 
the grids. 

2. We use a cutoff dipole interaction potential for im- 
proved accuracy and to reduce the size of the numerical 
grids needed. The cutoff potential minimizes interac- 
tion between aliased copies of the system (problematic 
with Fourier methods used for systems with long-range 
interactions). We make use of both the cylindrical cut- 
off suggested in [12] and the spherical version devel- 
oped in [17]. 

3. We use two grid extents as schematically shown in 
Fig. 7(b). Since we only calculate the discrete modes up 
to some relatively small energy, e cut , we can use a small 
and dense set of grids to accurately perform the diago- 
nalizations [and obtain rid(x)]. Then a much larger grid 
is used for the semiclassical region which extends out to 
much higher energy, as needed to accurately capture the 
thermal tails of the system. 



4. In application to mechanical stability, finding self- 
consistent solutions is challenging and care needs to be 
taken to ensure that metastable states are not lost pre- 
maturely and that a coarse grid does not disguise insta- 
bility. In using fixed point iteration effectively, we em- 
ploy an algorithm to efficiently increase or decrease the 
mixing speed A mix during the self-consistent process. 
We have observed that if A m ; x is too large early on in 
self-consistency iterations a metastable solution may be 
lost. Normally we start with A m i x = 0.3 (<~ 0.01 for 
biconcave density regions) and appropriately increase 
or decrease this during the search for a self-consistent 
solutions depending upon conditions. We also note that 
care needs to be taken to reliably detect mechanical in- 
stability collapse. We perform a number of tests to de- 
termine instability including detecting the development 
of density spikes and large gaps in the low energy spec- 
trum. We have confirmed that these are good signatures 
of the grid dependent numerical collapse discussed in 
Sec. IIIB1. 



3. Convergence tests of stability boundary 

In oblate geometries the self consistent calculations above 
T c become increasingly difficult to perform as A increases. 
For this reason we have not extended our calculations beyond 
A = 8. The origin of the difficulty is two fold: (i) the effective 
potential flattens considerably which increases the low energy 
density of states, meaning that a large number of modes exist 
below e cut . (ii) In oblate geometries the interaction strengths 
at instability are much larger, and these numerous low energy 
modes interact strongly with each other. 

The important convergence test is that our results are inde- 
pendent of the boundary (e cut ) separating the discrete low en- 
ergy modes from the continuous semiclassical spectrum. The 
typical error bars shown in Figs. 3 and 4 represent the varia- 
tion in the stability boundary as e cut was varied (1 a spread 
obtained from tests where the number of discrete modes 
ranged over <~ 10 to 1000). It is not computationally feasi- 
ble for us to take e cut high enough for A > 4 (above T c ) to 
get a self-consistent result fully independent of e cut . How- 
ever, note we do not observe a systematic drift nor monotonic 
relationship between e cut and the stability boundary. 

Below T c the error bars are instead determined arbitrarily 
by the bisection tolerance of parameters [i and dd- 
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